home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / LIBIP / N2_PV.C < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  11.5 KB  |  493 lines
  1. /* 
  2.  * n2_pv.c
  3.  * 
  4.  * Practical Algorithms for Image Analysis
  5.  * 
  6.  * Copyright (c) 1997, 1998, 1999 MLMSoftwareGroup, LLC
  7.  */
  8.  
  9. /*
  10.  * N2_P(oly)V(ert)
  11.  * Misc polygon/vertices functions
  12.  *
  13.  * find 2**n new vertices by resampling input polygon 
  14.  * 
  15.  * determine radius vector relative to centroid; shift to zero mean;
  16.  * reconstruct angular bend function with respect to new set of vertices;
  17.  *
  18.  */
  19.  
  20. #include <stdio.h>
  21. #include <stdlib.h>
  22. #include <math.h>
  23. #include "ip.h"
  24.  
  25. #define ZERO        0.0
  26.  
  27. #define    F_TO_INT(x)    ( ((x)-(int)(x)<0.5) ? (int)(x) : (int)(x)+1 )
  28. #define    SQR(a)        ( (a)*(a) )
  29. /*#define PERP(a, b)    ( -((a)->y)*(b)->x+(a)->x*(b)->y ) */
  30.  
  31.  
  32. #define LOOP_COMPLETE    8.0
  33.  
  34. #define ON        1
  35. #define OFF        0
  36. #define SHIFT_CENTROID    ON
  37. #define    INQ_SHIFT_CTR    OFF
  38. #define    X_SHIFT        0.4
  39. #define    Y_SHIFT        0.4
  40.  
  41.  
  42. #define    DRAW_NEW_PTS    ON
  43. #define    RAD_CIRC    3
  44.  
  45. #define    ZERO_MEAN    0
  46. #define    P_AREA        1               /* guarantees vanishing spectrum for circs */
  47. #define    R_NORM        P_AREA          /* method of normalization for r_dev() */
  48.  
  49. #define    CHECK_CIRC
  50. #undef    D_DEBUG
  51.  
  52.  
  53. /*
  54.  * tg_turn_an()
  55.  *   DESCRIPTION:
  56.  *     evaluate tangent of the angle subtended by segments with
  57.  *     respective slopes m1, m2; two segments are perpendicular if m1*m2 = -1.0;
  58.  *   ARGUMENTS:
  59.  *     m1: first slope (double)
  60.  *     m2: second slope (double)
  61.  *   RETURN VALUE:
  62.  *     tangent between m1 and m2 (double)
  63.  */
  64.  
  65. double
  66. tg_turn_an (double m1, double m2)
  67. {
  68.   double tg_th, den;
  69.  
  70.   den = m1 * m2 + 1.0;
  71.   if ((-1.0e-10 < den) && (den < 1.0e-10))
  72.     tg_th = 500.0;
  73.   else
  74.     tg_th = (m1 - m2) / den;
  75.  
  76.   return (tg_th);
  77. }
  78.  
  79. /*
  80.  * slope()
  81.  *   DESCRIPTION:
  82.  *     evaluate slope of linear segment
  83.  *   ARGUMENTS:
  84.  *     pt1: first point (spoint *)
  85.  *     pt2: second point (spoint *)
  86.  *   RETURN VALUE:
  87.  *     slope between pt1 and pt2 (double)
  88.  */
  89.  
  90. double
  91. slope (struct spoint *pt1, struct spoint *pt2)
  92. {
  93.   double m;
  94.   double max_slope = 500.0;     /* max. poss. slope */
  95.   int x1, y1, x2, y2;
  96.  
  97.   x1 = F_TO_INT (pt1->x);
  98.   y1 = F_TO_INT (pt1->y);
  99.   x2 = F_TO_INT (pt2->x);
  100.   y2 = F_TO_INT (pt2->y);
  101.  
  102.   if (F_TO_INT (x2 - x1) == 0) {
  103.     if (F_TO_INT (y2 - y1) == 0)
  104.       m = 0.0;
  105.     else
  106.       m = max_slope;
  107.   }
  108.   else
  109.     m = ((double) (y2 - y1)) / ((double) (x2 - x1));
  110.  
  111.   return (m);
  112. }
  113.  
  114. /*
  115.  * vert_to_phi_star()
  116.  *   DESCRIPTION:
  117.  *     compute cumulative angular bend function for set of vertices
  118.  *     by evaluating turn angles of contour at each vertex
  119.  *   ARGUMENTS:
  120.  *     v: vertices (spoint *)
  121.  *     nv: # vertices (long)
  122.  *     phi_star: turn angles (float *)
  123.  *   RETURN VALUE:
  124.  *     currently always 1
  125.  */
  126.  
  127. int
  128. vert_to_phi_star (struct spoint *v, long nv, float *phi_star)
  129. {
  130.   //int i;
  131.   //float sign;
  132.   float ang_step;
  133.   float eps = (float) 0.5;
  134.  
  135.   printf ("\n...new version of function not complete\n");
  136.   return (-1);
  137.  
  138.   ang_step = (float) (2.0 * PI / (float) nv);
  139.  
  140. /*
  141.  * compute contour turn angle at point (xi, yi) from the 
  142.  * triple of points (ximm, yimm), (xi, yi) and (xipp, yipp)
  143.  * and collect results in phi_k;
  144.  * sign convention:
  145.  *   CW turn is assigned a negative turn angle
  146.  *
  147.  * (see poly_app.c, hsb)
  148.  */
  149.  
  150.  
  151. /* 
  152.  * check for closure 
  153.  * 
  154.  * ang_offset = *(n_phi_k + 0);
  155.  * for (i=0; i<=nv; i++)
  156.  * *(n_phi_k + i) -= ang_offset;
  157.  * 
  158.  * cum_turn_angle = *(n_phi_k + nv);
  159.  * printf("\n...cum_turn_angle = %f\n", cum_turn_angle);
  160.  * if( cum_turn_angle > 0.0 ) sign = -1.0;
  161.  * if( cum_turn_angle < 0.0 ) sign = 1.0;
  162.  * if( fabs(cum_turn_angle ) - 2.0*PI > eps )
  163.  * printf ("\n...total ang. bend not equal to two pi!!\n");
  164.  * 
  165.  * printf ("...sign = %f, ang_step = %f, ang_offset = %f\n",
  166.  * sign, ang_step, ang_offset);
  167.  */
  168.  
  169. /*
  170.  * evaluate phi_star by properly normalizing phi_k
  171.  * 
  172.  * for(i=0; i<=nv; i++)
  173.  * *(phi_star+i) = *(n_phi_k+i) + sign*((float)i)*ang_step;
  174.  */
  175.   return (+1);
  176. }
  177.  
  178.  
  179. /*
  180.  * zero_mom()
  181.  *   DESCRIPTION:
  182.  *     evaluate polygon area from given set of vertices as moment of order zero
  183.  *   ARGUMENTS:
  184.  *     v: vertices (spoint *)
  185.  *     nv: # vertices (long)
  186.  *   RETURN VALUE:
  187.  *     polygon area (double)
  188.  */
  189.  
  190. double
  191. zero_mom (struct spoint *v, long nv)
  192. {
  193.   int i;
  194.   double m00;
  195.   double ximm, xi, yimm, yi;
  196.  
  197.   m00 = 0.0;
  198.   for (i = 1; i <= (int) nv; i++) {
  199.     ximm = (double) (v + i - 1)->x;
  200.     xi = (double) (v + i)->x;
  201.     yimm = (double) (v + i - 1)->y;
  202.     yi = (double) (v + i)->y;
  203.  
  204.     m00 += yi * ximm - xi * yimm;
  205.   }
  206.   return (0.5 * fabs (m00));
  207. }
  208.  
  209.  
  210. /*
  211.  * shape_parm()
  212.  *   DESCRIPTION:
  213.  *     evaluate global shape parameter
  214.  *   ARGUMENTS:
  215.  *     c_len: curvature length (double)
  216.  *     moments: moments array (float *)
  217.  *   RETURN VALUE:
  218.  *     shape parameter (double)
  219.  */
  220.  
  221. double
  222. shape_parm (double c_len, float *moments)
  223. {
  224.   return (c_len * c_len / (4.0 * PI * (*(moments + 1))));
  225. }
  226.  
  227.  
  228. /*
  229.  * E_dist()
  230.  *   DESCRIPTION:
  231.  *     evaluate Euclidean distance between two points
  232.  *   ARGUMENTS:
  233.  *     pt1: first point (spoint *)
  234.  *     pt2: second point (spoint *)
  235.  *   RETURN VALUE:
  236.  *     distance between pt1 and pt2 (double)
  237.  */
  238.  
  239. double
  240. E_dist (struct spoint *pt1, struct spoint *pt2)
  241. {
  242.   double delx, dely, d;
  243.  
  244.   if (fabs (delx = (double) (pt2->x - pt1->x)) < 1.0)
  245.     delx = 0.0;
  246.   if (fabs (dely = (double) (pt2->y - pt1->y)) < 1.0)
  247.     dely = 0.0;
  248.   d = sqrt (SQR (delx) + SQR (dely));
  249.  
  250. #ifdef D_DEBUG
  251.   printf ("...d( (%3d,%3d), (%3d,%3d) ): %lf\n",
  252.           pt1->x, pt1->y, pt2->x, pt2->y, d);
  253. #endif
  254.   return (d);
  255. }
  256.  
  257. /*
  258.  * E_dist()
  259.  *   DESCRIPTION:
  260.  *     evaluate vector (cross) product of two input vectors
  261.  *   ARGUMENTS:
  262.  *     pt1: first point (spoint *)
  263.  *     pt2: second point (spoint *)
  264.  *   RETURN VALUE:
  265.  *     cross product between pt1 and pt2 (double)
  266.  */
  267.  
  268. double
  269. vcp (struct spoint *v1, struct spoint *v2)
  270. {
  271.   return (-(v1->y) * (v2->x) + (v1->x) * (v2->y));
  272. }
  273.  
  274.  
  275. /*
  276.  * cont_bend_en()
  277.  *   DESCRIPTION:
  278.  *     determine bending energy of closed contour by summing the square
  279.  *     of the curvature, based on the evaluation of the local radius of 
  280.  *     curvature, R, as the radius of the osculating circle, defined by 
  281.  *     a triplet of successive points on the contour; for three such points, 
  282.  *     v0, v1, v2, the radius R is given as R = d(v2, v0)/(2*sin th), where
  283.  *     th denotes the angle subtended by the segments s(v1, v0) and s(v2, v1),
  284.  *     and it is related to the contour turn angle by PI-th;
  285.  *
  286.  *     ref: D. J. Struik, ``Lectures on Classical Diff. Geometry'', Ch. 1-4
  287.  *         Addison-Wesley Publ. Co, 1961;
  288.  *   ARGUMENTS:
  289.  *     v: array of points (spoint *)
  290.  *     n: length of array (long)
  291.  *   RETURN VALUE:
  292.  *     bending energy of the cloaed contour (double)
  293.  */
  294.  
  295. double
  296. cont_bend_en (struct spoint *v, long n)
  297. {
  298.   int i;
  299.   double sn_th, dm, dp;
  300.   double kap, e_bend;
  301.   struct spoint *pvmm, *pv, *pvpp;
  302.   struct spoint sm, *psm = &sm, sp, *psp = &sp;
  303.  
  304.   e_bend = 0.0;
  305.   for (i = 0; i < n; i++) {
  306.     pv = &v[i];
  307.     if (i == 0)
  308.       pvmm = &v[n - 1];
  309.     else
  310.       pvmm = &v[i - 1];
  311.     if (i == n - 1)
  312.       pvpp = &v[0];
  313.     else
  314.       pvpp = &v[i + 1];
  315.  
  316.     psm->x = pv->x - pvmm->x;
  317.     psm->y = pv->y - pvmm->y;
  318.     psp->x = pvpp->x - pv->x;
  319.     psp->y = pvpp->y - pv->y;
  320.  
  321. /* 
  322.  * evaluate square of local curvature
  323.  */
  324.     dm = E_dist (pv, pvmm);
  325.     dp = E_dist (pvpp, pv);
  326.     sn_th = vcp (psm, psp) / (dm * dp);
  327.     kap = 2.0 * (sn_th / E_dist (pvpp, pvmm));
  328.  
  329.     printf ("...kap(%3d): %f\n", i, (float) kap);
  330.  
  331.     e_bend += (double) (SQR (kap) * 0.5 * (dm + dp));
  332.   }
  333.   return (e_bend);
  334. }
  335.  
  336.  
  337.  
  338. /*
  339.  * r_dev()
  340.  *   DESCRIPTION:
  341.  *     evaluate the function (r(s) - rbar)/rbar for poly {v} with centroid vc;
  342.  *
  343.  *     note:
  344.  *      ---> the radial function can become multivalued if "overhangs" occur,
  345.  *      even if the contour is a Jordan curve (enclosing a singly
  346.  *      connected region)!!
  347.  *      ---> the centroid does not in general coincide with the point with
  348.  *      respect to which r(s) must be evaluated to ensure
  349.  *      the vanishing of the first Fourier component of (r(s) - rbar)!!
  350.  *   ARGUMENTS:
  351.  *     vc: centroid (spoint *)
  352.  *     v: polygon vertices (spoint *)
  353.  *     r: (r(s) - rbar)/rbar (float *)
  354.  *     nv: # vertices in v (long)
  355.  *     moments: (float *)
  356.  *     imgIO: output image (Image *)
  357.  *     value: grayscale for drawing 0-255 (int)
  358.  *   RETURN VALUE:
  359.  *     none
  360.  */
  361.  
  362. void
  363. r_dev (struct spoint *vc, struct spoint *v,
  364.        float *r, long nv, float *moments, Image * imgIO, int value)
  365. {
  366.   int i;
  367.   float xfc, yfc;
  368.   float xcc = (float) X_SHIFT, ycc = (float) Y_SHIFT;
  369.   float rbar = (float) 0.0;
  370.   double delx, dely;
  371.   double m00 = -1.0;
  372.  
  373.   xfc = (float) vc->x;
  374.   yfc = (float) vc->y;
  375.   printf ("\nR_DEV: centroid (%3d, %3d)\n", vc->x, vc->y);
  376. /* 
  377.  * shift centroid coordinates to minimize first Fourier mode 
  378.  * (determine required shift by trial and error)
  379.  */
  380.   if (SHIFT_CENTROID == ON) {
  381.     if (INQ_SHIFT_CTR == ON) {
  382.       printf ("\n...current centroid position (xc, yc)\n");
  383.       printf ("-->correct xc by:");
  384.       scanf ("%f", &xcc);
  385.       printf ("-->correct yc by:");
  386.       scanf ("%f", &ycc);
  387.     }
  388.     xfc += xcc;
  389.     yfc += ycc;
  390.     printf ("...corrected centroid position: (%f, %f)", xfc, yfc);
  391.   }
  392.  
  393.   rbar = (float) 0.0;
  394.   for (i = 0; i < nv; i++) {
  395.     delx = (v + i)->x - xfc;
  396.     dely = (v + i)->y - yfc;
  397.     *(r + i) = (float) sqrt (delx * delx + dely * dely);
  398.  
  399.     rbar += *(r + i);
  400.   }
  401.  
  402.  
  403.   switch (R_NORM) {
  404.   case ZERO_MEAN:              /* shift r to zero mean and normalize */
  405.     rbar /= (float) nv;
  406.     break;
  407.   case P_AREA:                 /* compare with circle of equal area */
  408.     rbar = (float) sqrt (*(moments + 1) / PI);
  409.     break;
  410.   }
  411.   for (i = 0; i < nv; i++)
  412.     *(r + i) = (float) ((*(r + i) / rbar) - 1.0);
  413.  
  414. #ifdef CHECK_CIRC
  415.   printf ("\n...reference circle (method: %d):", R_NORM);
  416.   printf (" rbar: %f\n", rbar);
  417.  
  418.   draw_circle ((int) xfc, (int) yfc, (int) rbar, imgIO, value);
  419. #endif
  420. }
  421.  
  422. /*
  423.  * n2_vert()
  424.  *   DESCRIPTION:
  425.  *     find 2**n new vertices by resampling input polygon 
  426.  *     construct new vertices (spaced in equi_distant steps along the contour)
  427.  *     by simple linear interpolation
  428.  *   ARGUMENTS:
  429.  *     v: input polygon (spoint *)
  430.  *     s: pointer to partial sums array (float *)
  431.  *     nv: # vertices in v (long)
  432.  *     v_n2: resampled polygon (spoint *)
  433.  *     nv2: length of v_n2 (long)
  434.  *     imgIO: output image (Image *)
  435.  *     value: grayscale for drawing 0-255 (int)
  436.  *   RETURN VALUE:
  437.  *     none
  438.  */
  439.  
  440. void
  441. n2_vert (struct spoint *v, float *s, long nv, struct spoint *v_n2, long nv2, Image * imgIO, int value)
  442. {
  443.   int i, i_n2;
  444.   double cl_in2, delta_cl;
  445.   double xi, ximm, delx;
  446.   double yi, yimm, dely;
  447.   double si, simm, dels;
  448.  
  449.   double eps = 0.01;
  450.  
  451.   int c_index = 127;
  452.  
  453.  
  454.   (v_n2 + 0)->x = (v + 0)->x;
  455.   (v_n2 + 0)->y = (v + 0)->y;
  456.   delta_cl = 100.0 / ((double) nv2);
  457.   for (i_n2 = 0; i_n2 <= nv2; i_n2++) {
  458.     i = 1;
  459.     cl_in2 = ((double) i_n2) * delta_cl;  /*norm. cont. len */
  460.     while ((cl_in2 - *(s + i)) > eps)
  461.       i++;                      /* find edge */
  462.  
  463.     /* new vertex may be close to an existing one */
  464.     if (fabs (*(s + i) - cl_in2) <= eps) {
  465.       (v_n2 + i_n2)->x = (v + i)->x;
  466.       (v_n2 + i_n2)->y = (v + i)->y;
  467.     }
  468.     /* in general, it won't be: interpolate */
  469.     else if (fabs (*(s + i) - cl_in2) > eps) {
  470.       xi = (double) (v + i)->x;
  471.       ximm = (double) (v + i - 1)->x;
  472.       delx = xi - ximm;
  473.       yi = (double) (v + i)->y;
  474.       yimm = (double) (v + i - 1)->y;
  475.       dely = yi - yimm;
  476.       si = *(s + i);
  477.       simm = *(s + i - 1);
  478.       dels = si - simm;
  479.  
  480.       (v_n2 + i_n2)->x = F_TO_INT (delx * (cl_in2 - simm) / dels + ximm);
  481.       (v_n2 + i_n2)->y = F_TO_INT (dely * (cl_in2 - simm) / dels + yimm);
  482.     }
  483.   }
  484.  
  485. /* 
  486.  * overlay new nv2 vertices on displayed polygon 
  487.  */
  488.   if (DRAW_NEW_PTS == ON) {
  489.     for (i_n2 = 0; i_n2 < nv2; i_n2++)
  490.       draw_circle ((int) (v_n2 + i_n2)->x, (int) (v_n2 + i_n2)->y, (int) RAD_CIRC, imgIO, value);
  491.   }
  492. }
  493.